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CONTRACTOR REPORT 


A TWO-LAYER MULTIPLE-TIME-SCALE TURBULENCE MODEL AND 

GRID INDEPENDENCE STUDY 

INTRODUCTION 


For wall bounded turbulent boundary layer flows, the logarithmic velocity pro- 
file prevails almost universally throughout the thin shear layer. The wall function 
method is based on the validity of the logarithmic velocity profile and the local 
equilibrium condition (the generation rate is almost equal to the dissipation rate of 
the turbulent kinetic energy) of turbulence. The wall function method lost its 
validity in many instances (i.e., separated flows, turbulent boundary layer flows 
subjected to strong adverse pressure gradient, and transition flows for which the 
local turbulent Reynolds number is sufficiently low). Nevertheless, partly due to 
the consideration of computational efficiency and partly due to the difficulty in 
developing and implementing the low Reynolds number turbulence models for elliptic 
flows, the wall function method has been used dominantly even for separated and 
recirculating flows. The appropriateness of the wall function method began to be 
re-evaluated [1], and a number of two-layer turbulence models [2-5] and a number 
of low Reynolds number turbulence models [6,7] began to appear in recent years. 

A comparative study on various low Reynolds number turbulence models can be found 
in Reference 7. 

In the simplest form of the two-layer turbulence models, the turbulent eddy 
viscosity inside the near- wall, low turbulent Reynolds number layer is expressed 
using the mixing length theory which is based on the Kolmogorov-Prandtl theory [8] 
of turbulence. A few flow cases in which the simplest form two-layer model yielded 
significantly improved computational results over the wall function method can be 
found in Reference 1. The two-layer turbulence model presented in this report 
belongs to the simplest class of two-layer models in a sense that the mixing length 
assumption has been used to derive the turbulent eddy viscosity expression for the 
near wall low turbulent Reynolds number region. 

In some of the finite element computation of turbulent flows [9-11], the tan- 
gential velocity has been evaluated using the wall function method and the vanishing 
normal velocity has been prescribed at the near wall boundary. Use of the vanishing 
normal velocity at the near wall boundary has been justified based on the assumption 
that the magnitude of the normal velocity would be vanishingly small. However, it 
brought about an uncertainty with regard to satisfying the conservation of mass 
equation, since the flow below the near wall boundary (the viscous sublayer and the 
transition region) may accelerate or decelerate depending on the pressure distribution 
in the region, and accordingly, the normal velocity component may seldom vanish. 
Slightly different implementations of the wall function method to satisfy the conserva- 
tion of mass constraint more rigorously can be found in References 12 and 13. In 
Reference 12, the normal velocity at the near wall boundary has been evaluated by 
integrating the conservation of mass equation in the region below the near wall 
boundary; and in Reference 13, the normal velocity at the near wall boundary has 
been estimated using an analytical expression obtained from a near-wall analysis based 
on mixing length theory with the pressure gradient effect excluded. Implementation 
of such methodologies for three-dimensional flows may be complicated or exceedingly 
difficult. In many finite difference methods, the difficulty of satisfying the conserva- 
tion of mass equation has been circumvented by using the staggered grids [ 14] . 



In the two-layer model presented herein, the computational domain for the flow 
equations penetrated up to the wall, where no-slip boundary condition has been 
prescribed. Use of the two-layer model in a finite element computation of turbulent 
flows partly eliminated an uncertainty involved in satisfying the conservation of mass 
constraint. Grid independence study has been made for the near-wall layer; and the 
computational results have been compared with the experimental data as well as those 
obtained by using the wall function method. The purpose of the present study has 
also been to clarify the uncertainty in satisfying the conservation of mass equation, 
which has been caused by using the vanishing normal velocity in the context of the 
wall function method, by comparing the computational results obtained by using the 
two-layer model with those obtianed by using the wall function method. 

In finite difference methods, it is known that a significant amount of grid 
points are required inside the near-wall layer. Discussion on computational efficiency, 
in comparison with the finite difference methods, is also included. 


NEAR WALL EDDY VISCOSITY EQUATION FOR TWO-LAYER MODEL 


Very close to the wall, inside the viscous sublayer, the eddy viscosity grows 
in proportion to the cubic power of distance from the wall and the mean tangential 
velocity varies in proportion to the same wall distance. The eddy viscosity equation 
given below has been derived by incorporating these experimentally observed behavior 
of turbulence into the logarithmic velocity profile equation, and the modified logarithmic 
velocity profile equation has been derived by integrating the eddy viscosity equation. 
Detailed derivation of these equations are given in Reference 15. 


■> t + = i E {exp( K u + > - [1 + (ku + ) + \ (ku + ) 2 ]} 


( 1 ) 
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where v t + (v t + = v t /v) is a non-dimensional turbulent eddy viscosity; v t is the turbu- 
lent eddy viscosity; v is the kinematic viscosity of the fluid; u + (u + = u/u^) is a 


non-dimensional velocity; u 
wall shearing stress; y + (y 
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[u = /(t w /p)] is the wall friction velocity; t w is the 
= u^y/v) is the wall coordinate; k is the von Karmann 


constant; and E is an experimentally determined constant coefficient . The values of 
k and E used in Reference 15 are given as 0.4 and 7.7, respectively; whereas in most 
of numerical analysis of turbulent flows using the k-e type turbulence models, k = 

4.1 and E = 9 have been used. The eddy viscosity and the modified logarithmic 
velocity profiles obtained by using these two different sets of coefficients are com- 
pared with experimental data in Figures 1 and 2. It can be seen that these two sets 
of coefficients yielded almost identical results. In the present study, k = 4.1 and 
E = 9 have been used throughout. 
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HIGH REYNOLDS NUMBER TURBULENT BOUNDARY LAYER FLOW EQUATIONS 


The turbulent boundary layer flow equations and the high turbulent Reynolds 
number multiple-time-scale turbulence model are described briefly below for complete- 
ness. Details on the turbulence model can be found in Reference 11. 

The turbulent boundary layer flow equations are given as: 
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where u and v are the time averaged mean velocities in flow direction and in trans- 
verse direction, respectively, and p is the pressure. The turbulent eddy viscosity 
inside the near-wall layer for the mean momentum equation has been computed by 
using equation (1). The turbulent eddy viscosity outside the near-wall layer has 
been obtained by using the eddy viscosity equation given as [ 11] , 



(5) 


where c - (=0.09) is a constant, k (k = k + k.) is the turbulent kinetic energy, 
kp is the turbulent kinetic energy of eddies in the production range, k t is the turbu- 
lent kinetic energy of eddies in the dissipation range, and £p is the transfer rate of 
turbulent kinetic energy from the production range to the dissipation range. 

Transport equations for the turbulent kinetic energies for the multiple- time - 
scale turbulence model are given as: 


u 


9k p 9k a 

— 2. + v 2- - — r(v 

3 x 3y 3 y uv 


kp 


9 k 

) T ^]=Pr-c t 


u 


3k^ 

Tx~ 


+ v 


3k t 

9y 


_3 . t + V t , 3k t, _ 

9y I(V a kt ) 3 y e p 


(6) 


(7) 


where Pr is the production rate of turbulent kinetic energy, e t is the dissipation 
rate of the turbulent kinetic energy, and and are the turbulent Prandtl 
numbers for the turbulent kinetic energy equations. 
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The convection- diffusion equations for the energy transfer rate and the dissi- 
pation rate are given as: 
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where c^ (£ = 1,3) and c^ (£ = 1,3) are the turbulence model constants, and o 

and a ^ are the turbulent Prandtl numbers for the energy transfer rate and the 

energy dissipation rate equations, respectively. The turbulence model constants are 
given as: = 0.75, = 1.15, = 0.75, = 1.15, Cp^ = 0.21, Cp 2 = 1.24, 

Cpg = 1.84, c^ = 0.29, c ^2 = 1.28, and c^g = 1.66 [11]. 

WALL FUNCTION BOUNDARY CONDITIONS 


The near wall boundary conditions for k p , k t , e p , and e t were obtained from 

the standard wall function method. Detailed derivation of the wall function boundary 
conditions for the multiple -time- scale turbulence model can be found in Reference 11. 
These are given as. 
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COMPUTATIONAL RESULTS 


The finite element computational procedure, as related to the two-layer approach, 
is described briefly below. 

The governing differential equations were solved on the physical domain using 
physical dimensions. The wall shearing stress to be used in equations (1), (11), and 
(12) has been obtained by using equation (10). Alternatively, the wall shearing stress 
could have been computed by using the standard definition given as, 
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3 yjy=0 


(15) 


Equation (15) is valid at and very close to the wall where the molecular viscosity 
dominates over the turbulent eddy viscosity. However, an accurate evaluation of the 
velocity gradient at the wall may require a significant number of grid points inside 
the near-wall layer. Furthermore, the turbulent kinetic energies, the energy transfer 
rate, and the dissipation rate were computed using the standard wall function method. 
Therefore, equation (10) has been preferred over equation (15) for the two-layer 
model presented herein. 

All the governing differential equations and the wall function boundary condi- 
tions are nonlinearly coupled with each other; thus, the system of equations, including 
the wall function boundary conditions, have been solved iteratively until the prescribed 
convergence criterion was satisfied. A direct (Picard) iteration method has been used 
to solve the system of equations [10]. Details on the finite element method for turbu- 
lent flows can be found in Reference 10. 

The initial condition data for velocity, turbulent kinetic energy, and dissipation 
rate were obtained from experimental data. The initial condition data for the ratios 
of kp/kj and £p/^ were obtained by interpolating the near wall values and the free 

stream values of these ratios, respectively [11]. The wall function boundary condi- 
tions, equations (10) through (14), were used at the near wall boundary for the tur- 
bulence equations; and the vanishing gradient boundary condition was used at the 
outer edge of the computational domain for both the flow equations and the turbulence 

equations. Whenever necessary, the values of 1.225 kg/m and 0.17854 x 10 kg/ 
m-sec were used for density and molecular viscosity, respectively. 


1. A Wall Jet Issuing Into a Moving Stream 

The configuration of the wall jet flow considered herein can be found in Refer- 
ences 10 and 18, and the experimental data for the flow can be found in Reference 18. 

Input data used in computation of the flow were obtained directly and/or by 
curve-fitting the experimental data [10,18]. The transverse computational domain of 
the near-wall layer extended from the wall (y = 0 m) to y = 0.002 m; and that of the 
outer layer extended from y = 0.002 m to y = 0.08 m. The computational domain in 
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the flow direction extended from x = 0.5532 m (x/b = 82.2, where b is the jet slot 
width) to x = 1.6892 m (x/b = 251). Two different discretizations have been used in 
the present study. In the first case, the near-wall layer has been discretized by 
using 20 grids (10 equally space quadratic elements); and in the second case, two 
grid points (1 quadratic element) have been allocated in the near-wall layer. For 
both of the cases, the transverse domain in the outer layer has been discretized by 
45 unequally spaced quadratic elements, and the flow direction domain has been dis- 
cretized by 1135 line-steps. For both of the cases, an average of 12 iterations was 
required for each line- step to obtain convergent solutions. This number of iterations 
was found to be the same as the one given in the previous study [11], in which the 
standard wall function method had also been used for the flow equations at the near 
wall boundary. 

The computational results of the flow development along the downstream direc- 
tion, obtained by using the two-layer model, are compared with experimental data as 
well as those obtained by using the standard wall funciton method in Figure 3(a). 

The velocity, turbulent kinetic energy, and Reynolds stress at the far downstream 
location (x/b = 251) are shown in Figures 3(b); and the computed mean normal 
velocity profiles at the same downstream location are shown in Figure 3(c). 

It can be seen in Figures 3(a) and 3(b) that the computational results for the 
tangential mean velocity and the turbulence quantities obtained by using the two- 
layer model are almost the same as those obtained by using the standard wall function 
method. However, the mean normal velocity obtained by using the two-layer model 
was slightly different from the one obtained by using the wall function method 
[Fig. 3(c)]. For the wall-jet flow, the characteristic mean tangential velocity was a 
few orders of magnitude greater than that of the mean normal velocity, hence, the 
mean tangential velocity and the turbulence quantities were not significantly influenced 
by the mean normal velocity. 

It was found that the number of elements used to cover the near- wall layer did 
not exhibit significant influence on the converged solutions for all of the turbulent 
flow variables, and that the grid independent solution has been obtained with as small 
as two grid points (1 quadratic element) in the near wall layer. 


2. Wake-Boundary Layer Interaction Flows 

The configuration and the experimental data of the weakly coupled wake-boundary 
layer interaction flow considered herein can be found in Reference 19. 

The transverse domain of the near- wall layer extended from the wall (y = 0 m) 
to y - 0.002 m; and that of the outer layer extended from y = 0.002 m to y = 0.12 m 
in the free stream region. The computational domain in the flow direction extended 
from x = 0.2 m (x/d = 20, where d is the diameter of the cylinder submerged in the 
boundary layer flow) to x = 0.86 m (x/d = 86). Two different discretizations have 
been used for the near- wall layer, the details of which are the same as in the previous 
wall- jet case. The transverse domain in the outer layer has been discretized by using 
45 unequally spaced quadratic elements , and the flow direction domain has been dis- 
cretized by using 825 line steps. Approximately 10 iterations were required for each 
line-step for all of the cases, including the case of using the wall function method. 

The velocity, turbulent kinetic energy, and Reynolds stress at the far down- 
stream location (x/d = 86) obtained by using the two-layer model are compared with 
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experimental data as well as those obtained by using the wall function method in 
Figures 4(a), and the computed mean normal velocity profiles at the same downstream 
location are shown in Figure 3(b). Again, the computational results for the tangential 
mean velocity and the turbulence quantities obtained by using the two-layer model 
were almost identical to those obtained by using the wall function method [Fig. 4(a). 
The mean normal velocity obtained by using the two- layer model was slightly different 
from the one obtained by using the wall function method [Fig. 4(b)]. As in the 
wall- jet case, the mean tangential velocity and the turbulence quantities were not 
significantly influenced by the mean normal velocity, and the grid independent solu- 
tions were obtained with as small as two grid points inside the near- wall layer. 


CONCLUSIONS AND DISCUSSION 


A two-layer multiple-time-scale turbulence model and grid independence study 
for turbulent boundary layer flows have been presented in this report. In the two- 
layer model, the computational domain for the conservation of mass equation and the 
mean momentum equation has been located at the wall where no slip boundary condi- 
tion could be specified. Therefore, the conservation of mass constraint could have 
been satisfied more rigorously with the two-layer model than with the standard wall 
function method. 

The mean normal velocity obtained by using the two-layer model was slightly 
different from the one obtained by using the standard wall function method. However, 
the tangential mean velocity and the turbulence quantities were not influenced sig- 
nificantly by the slight difference in the mean normal velocity. Hence, it was found 
that the performance of the turbulence model has not been masked by use of the wall 
function method for these example flow cases. 

In the present finite element computation of the turbulent boundary layer flows, 
the grid independent solution has been obtained with only two grid points (1 quadratic 
element) inside the near-wall layer. In finite difference computations, it is usually 
known that at least eight grid points are required inside the near-wall layer for the 
simplest class two-layer models [1], The required number of iterations to obtain 
convergent solutions was not increased by the use of the two-layer model compared 
with that of using the wall function method. To summarize, the two-layer model 
turned out to be more advantageous over the wall function method in finite element 
computation of turbulent boundary layer flows. 
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: Eq. (1) WITH K = 0.41 AND E = 9, 

: Eq. (1 ) WITH k = 0.40 AND E = 7.7 

O : EXPERIMENTAL DATA [KLEBANOFF) 

X : EXPERIMENTAL DATA [LAUFER] 

Figure 1. Effective turbulent eddy viscosity 



O : EXPERIMENTAL DATA (KLEBANOFF] 
X : EXPERIMENTAL DATA [LAUFER] 


Figure 2. Logarithmic velocity profiles 
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# , . 1 : COMPUT. MAX. EXCESS VELOCITY (u 0 ), 

EXTERNAL VELOCITY (u e ). WALL SHEARING STRESS (l w ). AND LENGTH SCALE 
(y 1/2). RESPECTIVELY; A, X.A.O: EXPT u 0 . Wo,. iVZ, RESPECTIVELY; 

A, 0: 20 GRIDS (10 ELEMENTS) AND 2 GRIDS (1 ELEMENT) INSIDE THE NEAR- 
WALL LAYER. RESPECTIVELY; C: WALL FUNCTION METHOD. 

Figure 3(a). Development of the wall-jet flow along the flow direction. 



. , : COMPUT. RESULTS FOR VELOCITY (u). TURBULENT 

KINETIC ENERGY (k). AND REYNOLDS STRESS ( uV ). RESPECTIVELY; 

X. O. A : EXPT DATA FOR u, k. -uV, RESPECTIVELY; A. B : 20 GRIDS ( 1 0 
ELEMENTS) AND 2 GRIDS (1 ELEMENT) INSIDE THE NEAR-WALL LAYER. 

RESPECTIVELY; C: WALL FUNCTION METHOD. 

Figure 3(b). Velocity, turbulent kinetic energy, and Reynolds stress of the 

wall- jet flow at x = 1.6882 m. 
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Figure 4(a) . 

the 


Figure 4(b). 
flow at 
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Velocity, turbulent kinetic energy, and Reynolds stress for 
wake-boundary layer interaction flow, at x = 0.86 m. 
[Notations are the same as in Figure 3(b).] 



Normal velocity profile for the wake-boundary layer interaction 
x = 0.86 m. [Notations are the same as in Figure 3(c).] 


☆UX GOVERNMENT PRINTING OFFICE: I9M - ftft4$4/MI29 


TECHNICAL REPORT STANDARD TITLE PAGE 


1. REPORT NO. 2. GOVERNMENT ACCESSION NO. 

NASA CR- 183653 


4. TITLE AND SUBTITLE 

A Two-Layer Multiple-Time-Scale Turbulence Model and 
Grid Independence Study 


7. AUTHOR (S) 

JS.-W. Kim* and C. -P. Chen** 


9. PERFORMING ORGANIZATION NAME AND ADDRESS 

George C. Marshall Space Flight Center 
Marshall Space Flight Center, Alabama 35812 


12. SPONSORING AGENCY NAME AND AOORESS 

National Aeronautics and Space Administration 
Washington , D . C . 20546 


3. RECIPIENT'S CATALOG NO. 


5. REPORT DATE 




6. PERFORMING ORGANIZATION CODE 

ES42 


8. PERFORMING ORGANIZATION REPORT tt 


10. WORK UNIT NO. 

11. CONTRACT OR GRANT NO. 

NAS8-35918 

13. TYPE OF REPORT & PERI00 COVERED 

Contractor Report 


14. SPONSORING AGENCY CODE 


ts. supplementary notes Prepared for the Fluid Dynamics Branch, Space Science Laboratory 

Technical Monitor : N. C. Costes 

♦Universities Space Research Association 

♦♦Department of Mechanical Engineering, UAH, Huntsville, Alabama 35899. 


16. ABSTRACT 

A two-layer multiple-time- scale turbulence model is presented in this report. 

The near-wall model is based on the classical Kolmogorov-Prandtl turbulence hypothesis 
and the semi- empirical logarithmic law of the wall. In the two-layer model presented 
herein, the computational domain of the conservation of mass equation and the mean 
momentum equation penetrated up to the wall, where no slip boundary condition has 
been prescribed; and the near wall boundary of the turbulence equations has been 
located at the fully turbulent region, yet very close to the wall, where the standard 
wall function method has been applied. Thus, the conservation of mass constraint can 
be satisfied more rigorously in the two-layer model than in the standard wall function 
method. In most of the two-layer turbulence models, the number of grid points to be 
used inside the near-wall layer posed the issue of computational efficiency. The present 
finite element computational results showed that the grid independent solutions were 
obtained with as small as two grid points, i.e., one quadratic element,^ inside the 
near- wall layer. Comparison of the computational results obtained by using the 
two-layer model and those obtained by using the wall function method is also presented. 


17. KEY WORDS 

Turbulence Model, Computational, Fluid, 
High Reynolds Number Flow, Finite 
Element 


18. DISTRIBUTION STATEMENT 


Unclassified — Unlimited 


E. Tandberg-faanssen 


19* SECURITY CLASSIF. (of thU report) 

20. SECURITY CLASSIF. (of thi» page) 

21. NO. OF PAGES 

22. 

PRICE 

Unclassified 

Unclassified 

17 


NTIS 


MSFC- Form 3393 (Rev. December It 73) 


For sale by National Technical Information Service, Springfield, Virginia 12161 

















